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ABSTRACT-  A hybrid  scheme  consisting  of  a modified 
second  order  in  time-  fourth  order  in  space  finite-difference 
time-domain  (FDTD)  scheme  ” M3d24 " and  the  Yee  algorithm, 
with  subgridding  is  introduced  to  overcome  the  errors  of 
applying  the  4th  order  in  space  FDTD  at  the  interfaces  of 
perfect  electric  conductors  (PEC)  or  dielectric  scatterers.  This 
hybrid  scheme  is  based  on  applying  the  Yee  algorithm  in  the 
vicinity  of  the  scatterer  using  a high  resolution  grid  (number  of 
points  per  wavelength),  and  the  M3d24  scheme  in  the  other 
regions  using  a low  resolution  grid  in  order  to  reduce  the 
required  computer  storage  for  large  problems,  while  still  good 
accuracy.  The  results  of  this  hybrid  scheme  are  shown  to  agree 
well  with  the  results  of  the  Yee  algorithm  using  a high 
resolution  grid,  for  problems  of  plane  wave  scattering  from 
PEC  cubes,  spheres. 
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subgridding,  transient  scattering  of  electromagnetic  waves. 

1.  INTRODUCTION 


storage  for  large  problems.  The  high  resolution  is  also  required 
to  reproduce  fast  variations  of  the  fields  in  certain  regions.  In 
addition,  complex  structures  which  contain  fine  details  in 
some  regions  need  high  resolution  grids,  whereas  a 4th  order 
scheme  with  low  resolution  can  be  applied  in  the  remaining 
regions.  A number  of  examples  are  given  in  the  present  paper 
to  show  the  applicability  of  the  scheme  for  different  categories 
of  problems. 


2.  METHOD  OF  SOLUTION 


The  modified  3D  4th  order  FDTD  scheme  " M3d24"  is 
introduced  by  the  authors  in  [2,3,4].  The  updating  equations  of 
the  M3d24  scheme  have  been  derived  as  applications  of  a 
version  of  Ampere’s  law  and  Faraday’s  law  using  2nd  order 
central  difference  in  time  and  4th  order  in  space. 

The  updating  equation  of  the  filed  component  Ex  is, 
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Higher  order  FDTD  algorithms  are  used  to  improve  the 
accuracy  of  the  computations  when  using  a relatively  low 
resolution  grid.  Taflove  [1]  approximated  the  derivatives  in  the 
conventional  fourth  order  FDTD  using  the  second  order 
central  difference  in  time  and  fourth  order  in  space.  A 
modified  second  order  in  time-  fourth  order  in  space  FDTD 
scheme  “A/3^24”,  has  been  deduced  by  the  authors  in  [2,3,4]. 
The  algorithm  enables  the  numerical  phase  velocity  error  to  be 
minimized,  so  that  it  leads  to  high  accuracy  with  low 
resolution  grids.  The  application  of  4th  order  algorithms 
directly  in  the  vicinity  of  a PEC  or  a dielectric  scatterer 
introduces  errors  due  to  using  large  stencils  in  the  4th  order 
algorithms  [5,6].  Other  sources  of  error  are  the  large  variations 
of  the  fields  in  the  vicinity  of  an  edge  of  a PEC  or  a large 
curvature  of  a surface  of  a PEC,  that  can  not  be  reproduced 
with  low  resolution.  To  overcome  the  error  of  large  stencils, 
Hadi  [5]  applied  the  2nd  order  Yee  algorithm  in  the  vicinity  of 
the  scatterer,  however,  the  Yee  algorithm  produces  a 
significant  error  if  it  is  applied  using  low  resolution. 
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The  approach  adopted  in  this  paper  depends  on  the 
application  of  the  Yee  algorithm  in  the  vicinity  of  the  scatterer 
using  a high  resolution  grid,  and  using  theA/3d24  scheme  with 
low  resolution  in  the  other  regions  to  reduce  the  required 
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where  each  point  in  the  Yee  cell  P(x,,y zt)  = P{ih,jh,kh ) 
and  the  time  tn  — nAt , h and  A?  are  the  space  and  time 

increments,  i,j,k,tl  are  integers  and  /y  = — ■ 

' eBh 

where  the  following  conditions  must  be  satisfied. 

Kl+K2+2(K}+KA)  = l (2) 

K,  K3_l-  (Kl  + Kl)_  (3) 

K2  AT,  2 (AT,  + /f2) 

The  constants  A|,A2,Al3and  are  optimized  to  minimize 
the  dispersion  error  in  all  directions  at  a certain  frequency  [3]. 
The  time  step  At  should  be  small  w.r.t  the  cell  size  according 
to  the  stability  condition  [3], 

At  < 7 (4) 

The  other  five  E,H  field  components  have  a similar  form 
corresponding  to  Ampere’s  and  Faraday’s  laws  in  the  different 
directions  with  the  field  locations  conforming  to  the  Yee  cell. 
The  magnetic  field  components  are  updated  at  times  shifted  by 
At  / 2 w.r.t  the  times  of  the  electric  field  components. 

The  present  approach  is  based  on  applying  the  M?>d24 
scheme  using  a coarse  grid  in  the  regions  which  are  not  in  the 
vicinity  of  objects,  and  applying  the  Yee  algorithm  using  a 
fine  grid  in  a region  that  contains  the  object  or  the  objects.  The 
ratio  of  the  resolution  of  the  fine  grid  to  that  of  the  coarse  grid 
to  be  used  equals  odd  integers.  An  important  property  of  using 
an  odd  integer  ratio  is  that  it  provides  colocated  fields  in  space 
for  the  fine  grid  with  the  coarse  grid  (see  Figure  1)  [7].  The 
time  step  used  is  entirely  constant  in  the  domain  with 

At  <A f , v (where  At  f is  the  maximum  allowable  time 

J max  J max 

step  satisfying  the  stability  condition  in  the  fine  grid). 

At  each  time  step  the  fields  are  updated  in  the  coarse  grid 
using  M3d24,  then  the  tangential  E-fields  and  the  normal  H- 
fields  on  the  interface  planes  of  the  fine  grid  are  interpolated 
from  the  coarse  grid  to  the  fine  one.  Finally  the  fields  are 
updated  in  the  fine  grid  using  the  Yee  algorithm. 


common  point  on  the  two  grids  with  indices  (7,y,  k)  on  the 
coarse  grid  and  on  the  fine  grid,  the  field 

Ev(i,j\k)  corresponds  to  the  field  ex  (/’+1,  j\  £’)  (noting 
that  the  position  of  Ex(i,jyk)  is  the  point 
((/  + 1 / 2)Ax, y*Ay, kAz\  Figure  (la)).  The  other  five  field 
components,  Ey(i,j,k)  corresponds  to  ey(i\j ’+1,^’)  , 
Ez(i,j,k)  corresponds  to  ez(i\j>k’+ 1)  , Hx(i,j\k ) 
corresponds  to  hx  {i\j ’+1,  k’+l)  , Hy(i,j,k ) corresponds 
to  hy(V+\,j\ k'+Y)  and  Hz(i,j,k)  corresponds  to 

hz(?+\,f+\,k'). 

In  order  to  show  the  actual  steps  of  evaluations  and 
interpolations  of  the  fields  in  the  coarse  and  the  fine  grids,  we 

consider  the  x-component  of  the  electric-field  Ex  and  ex. 

Assume  the  six  boundary  planes  of  the  fine  grid  region 
(intersection  planes  between  the  fine  and  coarse  grids)  are 

+i  = *'i4  Ax,  x2  = i2b  Ax , yl=jlbAy,  y2=j2bAy, 
Z|  = kib  Az  and  z2  — k2b  A z with  i2b  >7'u!  j2b  > J\b  and 
k2b  >klb . The  steps  of  the  computations  will  be  as  follows 
for  each  time  step; 

(1)  Ex  is  updated  in  the  coarse  grid  and  in  the  intersection 
planes,  i.e.  in  the  region  {(i2b  ^ i < ilb)  or 

( Jib  -J~  Jib  ) or  ( kib  ^ k ^ k\b  )}>  usinS  AT3 d2i  - 

(2)  Ex  is  interpolated  on  the  four  planes,  the  xy  planes 

Zj  — klh  Az  and  z2  = k2b  Az  and  the  xz  planes 

y\  = J\b  A y anc*  Y2  = J 2b  A y t0  ex  on  these 
planes.  The  interpolations  on  the  xy  plane 

Zj  = k lb  Az  (Figure  1),  can  be  evaluated  as  follows; 

i) Linear  interpolation  is  used  to  find  the  interpolated  value 
of  Ex  at  the  points  (p)  in  Figure  la,  which  are  the 

positions  of  the  ex  field  at  the  same  points  (i,j,k)  of  Ex . 

ii) The  interpolated  four  field  components  ex] ,ex4 

are  used  to  obtain  the  fields  at  the  points  (x)  (Figure  lb) 
using  quintic  interpolation. 

(3)  The  ex  field  is  updated  in  the  fine  grid  using  the  Yee 
algorithm. 

(4)  After  updating  the  ex  field,  the  values  of  Ex  are  updated 

as  E — ex  at  the  common  points  of  the  low  resolution 

and  the  high  resolution  grids  in  the  domain  of  the  fine 
grid. 

(5)  Simililar  steps  are  to  be  repeated  for  the  H field 
components. 


Assume  that  the  electric  and  magnetic  fields  are  denoted 
by  ( E , H ) in  the  coarse  grid  and  ( e , h ) in  the  fine  grid.  At  a 


The  use  of  the  Yee  algorithm  in  the  fine  grid,  rather  than 
using  the  M3d24  algorithm,  has  the  obvious  advantage  that  the 
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Yee  algorithm  does  not  have  a large  stencil  thus  it  introduces 
only  small  errors  at  the  object  boundary.  In  addition,  the  use  of 
high  resolution  will  increase  the  accuracy  of  the  Yee  algorithm 
to  the  required  value  while  using  fewer  floating  point 
operations  than  the  M3d1A  algorithm.  Third,  the  time  step 

A/ /max  ^ee  alS°rithm  *s  l°n&er  than  corresponding 

time  step  for  M3d2A  (according  to  the  stability  condition  of 
the  two  schemes  [1],  Eq.(4)),  i.e.  using  M3d24\r\  the  fine  grid 
will  increase  the  required  time  steps  to  solve  the  problem. 

3.  REFLECTION  ERROR  OF  SUBGRIDDING 

Using  multi-resolution  grids  in  the  computational  domain 
introduces  spurious  reflections  at  the  interfaces  between  the 
coarse  grid  and  the  fine  one.  These  reflections  are  evaluated 
using  a rectangular  waveguide  as  discussed  in  [7],  but  our  test 
is  applied  for  the  three-dimensional  case.  A TEl0  mode, 

Gaussian  pulse-modulated  sine  wave  is  excited  at  one  end  of  a 
waveguide  with  a cross  section  of  dimensions  20x9  mm  and 

its  center  lies  at  the  point  (*0,  J>0)  . The  total  field  is  updated 

taking  into  consideration  the  boundary  condition  of  the  perfect 
conductor  of  the  waveguide.  The  coarse  grid  has  a cell  edge 
length  h = Ax  = Ay  = tsz  =lmm.  The  waveguide  is  long 
enough  to  eliminate  reflections  from  the  absorbing  boundary 
condition  (ABC).  The  modulated  pulse  is  given  by 

g(t)=  sin  2 7if0t  e<t,T)l , (5) 


tY 


Figure  1.  Steps  of  interpolations  to  find  ex  from  E .in  the 
planes  z=£16Azor  z=k2b&z-  (a)  Interpolation  of  the  field 
values  at  the  points  (p),  (b)  the  interpolated  values  exl,ex2,ex3 
and  e are  used  to  obtain  the  field  values  at  the  points  (x). 

/0=1.5/Ciu=11.25GHZ, 

— = 20  GHZ  (6) 

It 

where  f q is  the  cuttoff  frequency  of  the  TEl0  mode. 

A fine  grid  region,  with  dimensions  16x4x40  mm  and  cell 
edge  length  h = Ax  = Ay  = Az  =l/3mm,  is  placed  such  that  its 
center  lies  in  the  transverse  plane  at  the  point  (x0 , y0  ) . It  is 

to  be  noted  that  the  reflections  are  due  to  all  six  interfaces  of 
the  subgridded  region  similar  to  the  test  region  of  Okoniewski 
et  al  [8],  which  must  be  greater  than  the  reflections  from  the  2- 
dimensional  test  structure  of  Chevalier  et  al  [7]  in  which  the 
fine  grid  filled  the  entire  second  half  of  the  waveguide. 

The  electric  field  component  £^is  computed  in  front 

of  the  subgridded  region  and  time  steps  are  taken  such  that  all 
reflections  from  the  subgridded  boundaries  are  taken  into 
consideration.  These  time  steps  are  stopped  before  any 
reflections  return  from  the  ABC.  The  same  computations  are 
repeated  using  entirely  the  coarse  grid,  i.e.  without 
subgridding,  and  the  resulting  field  values  are  taken  as 
reference  values.  The  difference  in  E between  the  two  cases 

represents  the  reflected  field  from  the  interfaces  of  the 
subgridded  region.  Figure  2a  shows  the  reflection  coefficient 
in  front  of  the  subgridded  region  at  the  centeral  point 

(*0,jk0)in  the  transverse  plane.  Figure2b  shows  the 

reflection  coefficient  as  an  average  of  the  reflection 
coefficients  for  different  x-positions  in  front  of  the  subgridded 
region  in  the  transverse  plane.  The  reflection  coefficients  in 
Figures  2a  and  2b  are  plotted  versus  frequency  from  7GHZ  to 


26 


ACES  JOURNAL,  VOL.  17,  NO.  1,  MARCH  2002,  SI:  APPROACHES  TO  BETTER  ACCURACY/RESOLUTION  IN  CEM 


43GHZ  (for  resolution  R=7  in  the  coarse  grid).  The 
subgridded  zone  is  expected  to  yield  multiple  reflections  along 
its  length  (40mm)  as  a half  wavelength  (corresponding  to 
3.75GHZ).  Peaks  at  approximately  multiples  of  this  frequency 
are  expected  to  occur,  Figures  2a  and  2b.  At  the  lower  end  of 
the  frequency  the  resolution  is  high  and  the  waveguide 
approaches  its  cut  off  frequency  which  leads  to  lower 
reflections.  It  is  seen  from  Figure  2 that  the  reflection 
coefficient  in  this  frequency  band  lies  below  -35dB  and  for 
resolution  R>10  ( f <30 GHZ)  the  reflection  coefficient  is 
smaller  than  -47dB. 


10.00  20.00  30.00  40.00  50.00 

f(GHZ) 


(a) 


10.00  20.00  30.00  40.00  50.00 


(b) 


f(GHZ) 


Figure  2.  Reflection  coefficient  from  a subgridded  region  16  x 
4 x 40  mm  in  a recatangular  waveguide  with  coarse  grid 
Ax=  Ay- Az  = l mm  and  refinement  factor  3.  (a)  Reflection 
coefficient  at  the  centeral  point  of  the  transverse  plane  in  front 
of  the  subgridding  region,  (b)  average  of  the  reflection 
coefficients  for  different  x-positions  in  the  transverse  plane. 


4.  APPLICATION  OF  THE  HYBRID 
SCHEME  TO  THE  FAR  FIELD 
CALCULATION  FROM  PLANE  WAVE 
SCATTERING  FROM  PEC  OBJECTS 

The  hybrid  scheme  is  tested  for  the  problem  of  scattering  from 
PEC  cubes  and  spheres  to  long  distances  to  show  the 
improvement  of  the  accuracy  compared  with  the  Yee 
algorithm.  Hadi  [6,5]  applied  the  2nd  order  Yee  algorithm  in 
the  vicinity  of  the  PEC  using  a single  resolution  in  the  entire 
domain.  The  application  of  the  Yee  algorithm  with  low 
resolution,  such  as  R=5,  introduces  a significant  error[5]. 

The  scattering  from  a cube  is  a good  test  to  show  the 
applicability  of  the  scheme  in  the  vicinity  of  PEC  scatterers 
which  contain  edges.  The  scattering  to  a long  distance  shows 
the  low  dispersion  error  of  the  scheme.  The  backscattered 
electric  field  from  a cube,  with  side  length  of  2 A A. , is  plotted 
in  Figure  3 with  the  distance  from  the  cube.  The  low 
resolution  (R=5)  results  of  various  schemes  are  compared  with 
the  high  resolution  (R=20)  result  of  Yee  algorithm.  It  is  seen 
from  Figure  3a  that  the  application  of  the  Yee  algorithm  with 
low  resolution  R=5  in  the  vicinity  of  the  PEC  cube  introduces 
a large  phase  error  which  persists  with  distance  and  affects  the 
solution  of  the  M3d2A  scheme  in  the  remaining  region.  In 
Figure  3b  the  use  of  the  M3d24  scheme  in  the  whole  region 
introduces  errors  in  the  amplitude  and  the  phase  due  to  using 
large  stencils  that  intersect  the  PEC  and  using  low  resolution 
near  the  edges  of  the  PEC.  The  result  of  applying  the  hybrid 
scheme,  in  which  the  Yee  algorithm  is  applied  with  R~15  in 
the  first  two  layers  around  the  cube  and  the  M3d2A  scheme  is 
applied  with  R=5  in  the  remaining  region,  is  shown  in  Figure 
3c.  It  is  seen  that  this  hybrid  scheme  gives  good  accuracy  in 
both  the  phase  and  the  amplitude  to  a long  distance. 

The  backscattered  electric  field  versus  distance  from 
PEC  spheres  is  shown  in  Figures  4 and  5,  for  spheres  with 
diameters  D = 1.2Aand  Z)  = 2.1/1,  respectively.  The  use  of 
low  resolution  near  the  PEC  sphere  introduces  a staircase  error 
which  can  be  reduced  by  using  high  resolution.  It  is  seen  from 
Figure  4a  that  the  result  of  Yee  algorithm  with  low  resolution 
R=7  does  not  give  accurate  results  in  the  vicinity  of  the  sphere, 
and  in  addition  a dispersion  error  appears  with  distance.  In 
Figure  4b,  the  subgridding  is  used  with  the  Yee  algorithm  such 
that,  on  the  surface  of  the  sphere  and  in  two  cells  around  the 
sphere,  the  Yee  algorithm  is  applied  with  resolution  R=21,  and 
the  same  algorithm  is  applied  with  R=7  in  the  remaining 
region.  The  application  of  the  Yee  algorithm  with  high 
resolution  in  the  vicinity  of  the  sphere  improves  the  results  in 
this  region,  and  the  low  resolution  introduced  a dispersion 
error  in  the  remaining  region,  as  seen  in  Figure  4b.  This 
dispersion  error  is  avoided  by  applying  M3d2A  without 
increasing  the  resolution  in  that  region  to  obtain  the  accurate 
results  of  Figure  4c.  Similarly,  it  is  seen  from  Figure  5 that  the 
hybrid  technique  M3c/24-Yee  with  R— 7/2 1 (the  Yee  algorithm 
is  applied  on  the  surface  of  the  sphere  and  in  two  cells  around 
the  sphere),  gives  accurate  results  in  amplitude  and  phase  to 
long  distance  from  the  sphere. 
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Figure  3.  The  backscattered  electric  field  with  distance  (dist  is 
the  distance  from  the  surface  of  the  cube)  for  normal  plane 
wave  incidence  on  a cube  (side  length 

Lx  =Ly  =4  = 2.4A  ),E,nc  =1000sin2^0r  jc,/0  =900MHZ, 
Yee  results  (R=20). 

(a) hybrid  M3d24-YEE  with  single  resolution  R=5,  (b) 

M3d24  with  R=5  (dotted  results),  (c)  hybrid  M3d24-YEE 
with  sugridding  R=5/15  (dotted  results) 


(a) 


(b) 


It  is  seen  from  Figures  3,4  and  5 that  the  electric  field 
has  reached  the  far  region  with  high  accuracy  when  using  the 
M3d24-YtQ  algorithm.  Thereby  this  technique  can  be  used  to 
find  directly  the  far-field,  with  a relatively  small  amount  of 
memory,  instead  of  using  a near  to  far  field  transformation. 
The  near  to  far  field  transformation  [9]  or  [10]  depends  on 
calculating  of  the  time  derivative  of  the  equivalent  electric  and 
magnetic  currents  which  affect  on  the  accuracy  as  shown  in 
[10].  An  improvement  is  achieved  in  [10]  by  implicitly 
calculating  the  required  time  derivatives  of  the  equivalent 
electric  and  magnetic  currents,  instead  of  an  explicit 
formulation  of  central  differencing  as  in  [9]. 

5.  CONCLUSION 

The  direct  application  of  4th  order  FDTD  algorithms  in  the 
vicinity  of  a perfect  electric  conductor  (PEC)  or  dilectric 
scatterer  introduces  errors  due  to  using  large  stencils  in  the  4th 
order  algorithms.  Other  sources  of  error  are  the  large 
variations  of  the  fields  in  the  vicinity  of  edges  of  a PEC  or  a 
large  curvature  of  a surface  of  a PEC  that  can’t  be  reproduced 
with  a low  resolution  in  the  entire  computational  domain 
because  of  a large  phase  error  due  to  the  application  of  the  Yee 
algorithm  with  low  resolution.  A hybrid  M3d24  - Yee  scheme 
with  subgridding  has  been  used  in  the  vicinity  of  the  scatterer 
using  a high-resolution  grid  and  M3d24  has  been  used  with  a 
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Figure  4.  The  scattered  electric  field  Esxcat  , with  distance,  from 
a sphere  with  diameter  1.2A.  (x0>y0)zQ)  is  the  center  of  the 

sphere,  E ""  = 1000sin2 jtf0tx,  f0  = 900 MHZ . 


Figure  5.  The  scattered  electric  field  with  distance  from  the 
surface  of  a sphere  with  diameter  2.1/1, 2?™  = 1000 sin  2nf0t  jc, 
f0  = 900MHZ  . 
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low-resolution  grid  in  the  other  region.  This  scheme 
successfully  overcomes  the  errors. 

Using  multi-resolution  in  the  computational  domain 
introduces  a spurious  reflection  at  the  interface  between  the 
coarse  grid  and  the  fine  one.  A ratio  of  three  has  been  used 
between  the  high-resolution  grid  and  the  low-resolution  grid. 
The  reflection  coefficient  of  the  subgridding  of  the  present 
algorithm  has  been  found  to  be  below  47  dB  for  frequences 
corresponding  to  resolution  greater  than  ten  in  the  coarse  grid. 

Different  test  cases  have  been  considered  which 
showed  the  applicability  of  this  hybrid  scheme  for  different 
types  of  problems  with  low  resolution.  The  scheme  has  been 
tested  for  scattering  from  PEC  cubes  and  spheres.  The 
scattering  from  a cube  is  a good  test  to  show  the  applicability 
of  the  scheme  in  the  vicinity  of  a PEC  scatterer  which  contains 
edges.  The  field  variation  with  distance,  up  to  the  far  field 
region,  shows  the  low  dispersion  error  of  the  scheme.  It  has 
been  found  that  the  field  reaches  the  far  region  for  the  test 
cases,  with  high  accuracy,  using  as  low  resolution  as  R=5  with 
the  present  scheme,  which  on  the  other  hand  requires  at  least  a 
resolution  R=20  in  the  Yee  algorithm.  Such  results  suggest  the 
use  of  this  scheme  to  obtain  the  far  field  directly  without  the 
need  to  use  near  field  to  far  field  transformation. 
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